pacman::p_load(sjPlot, jtools, dplyr, ggplot2, ggpubr, broom, pollster, tidyr, tidyverse, ggthemes, haven)

getwd()

gss <- read.csv("gss.csv")

# data cleaning
gss$pid <- car::recode(gss$partyid, '-99=NA; -98=NA; 0="Democrat"; 1="Democrat"; 2="Democrat"; 3="Independent"; 4="Republican"; 5="Republican"; 6="Republican"; 7=NA')
gss$partyid_without_leaners_as_partisans <- car::recode(gss$partyid, '-99=NA; -98=NA; 0="Democrat"; 1="Democrat"; 2="Independent"; 3="Independent"; 4="Independent"; 5="Republican"; 6="Republican"; 7=NA')

gss$republican <- ifelse(gss$pid=="Republican", 1, 0)
gss$democrat <- ifelse(gss$pid=="Democrat", 1, 0)
gss$independent <- ifelse(gss$pid=="Independent", 1, 0)

gss$republican_noleaners <- ifelse(gss$partyid_without_leaners_as_partisans=="Republican", 1, 0)
gss$democrat_noleaners <- ifelse(gss$partyid_without_leaners_as_partisans=="Democrat", 1, 0)
gss$independent_noleaners <- ifelse(gss$partyid_without_leaners_as_partisans=="Independent", 1, 0)

gss$age <- car::recode(gss$age, '-100=NA; -98=NA; -99=NA')
gss$age_cat_4 <- ifelse(gss$age<24.5, "18 to 24",
                        ifelse(gss$age>24.5 & gss$age<44.5, "25 to 44",
                               ifelse(gss$age>44.5 & gss$age<64.5, "45 to 64", "65 and over")))

gss$sex <- car::recode(gss$sex, '1="Male"; 2="Female"; -100=NA; -99=NA; -98=NA; -97=NA')
gss$sex <- factor(gss$sex, levels = c("Female", "Male"))

gss$race <- car::recode(gss$race, '-100=NA; 1="White"; 2="Black"; 3="Other"')

gss$region <- car::recode(gss$region, '1="New England"; 2="Middle Atlantic"; 3="East North Central"; 4="West North Central"; 5="South Atlantic"; 6="East South Central"; 7="Region_West South Central"; 8="Mountain"; 9="Pacific"')

gss$education <- car::recode(gss$degree, '-99=NA; -98=NA; -97=NA; 0="No College Degree"; 1="No College Degree"; 2="No College Degree"; 3="College Degree"; 4="College Degree"')
gss$degree <- car::recode(gss$degree, '-99=NA; -98=NA; -97=NA; 0="Less than High School"; 1="High School Degree"; 2="College Degree"; 3="College Degree"; 4="Graduate Degree"')
gss$degree <- factor(gss$degree, levels = c("Less than High School", "High School Degree", "College Degree", "Graduate Degree"))

gss$religion <- car::recode(gss$relig, '1="Protestant"; 2="Catholic"; 3="Other"; 4="Atheist/Agnostic/No Religion"; 5="Other"; 6="Other"; 7="Other"; 8="Other"; 9="Other"; 10="Other"; 11="Other"; 12="Other"; 13="Other"; -99=NA; -80=NA; -97=NA; -98=NA')
gss$religion <- factor(gss$religion, levels = c("Protestant", "Other", "Catholic", "Atheist/Agnostic/No Religion"))

gss$rur_urb <- car::recode(gss$xnorcsiz, '1="Urban"; 2="Urban"; 3="Suburban"; 4="Suburban"; 5="Suburban"; 6="Suburban"; 7="Rural"; 8="Rural"; 9="Rural"; 10="Rural"; -99=NA; -100=NA; -80=NA')
# City: 1,2; Suburb: 3,4,5,6; Rural: 9, 10; Small town: 7,8

gss$consci <- car::recode(gss$consci, '-100=NA; -98=NA; -99=NA; -97=NA; 1=3; 3=1')
gss$consci_binary <- ifelse(gss$consci==3,1,0)

gss$consci_cat <- car::recode(gss$consci, '1="Hardly any"; 2="Only some"; 3="A great deal"')
gss$consci_cat <- factor(gss$consci_cat, levels = c("Hardly any", "Only some", "A great deal"))
gss$consci_cat_binary <- ifelse(gss$consci_cat=="A great deal", "Percent 'A great deal' of trust in scientists", 0)
gss$consci_cat_lowest <- ifelse(gss$consci_cat=="Hardly any", "Percent 'hardly any' of trust in scientists", 0)

gss$attend <- car::recode(gss$attend, '-100=NA; -98=NA; -99=NA; -97=NA; 8="More than once a week";
                          7="Once a week"; 6="Once a week"; 5="A few times a month";
                          4="A few times a month"; 3="A few times a year"; 2="A few times a year";
                          1="Once a year or less"; 0="Never"')
gss$attend <- factor(gss$attend, levels = c("Never", "Once a year or less", "A few times a year",
                                            "A few times a month", "Once a week", "More than once a week"))

gss$class3 <- ifelse(gss$class_==1|gss$class_==2, "Lower/Working",
                     ifelse(gss$class_==3, "Middle",
                            ifelse(gss$class_==4, "Upper", NA)))

gss$class2 <- ifelse(gss$class_==1|gss$class_==2, "Lower/Working",
                     ifelse(gss$class_==3|gss$class3=="Upper", "Middle/Upper", NA))

gss$race_binary <- ifelse(gss$race=="White", "White",
                          ifelse(gss$race=="Black", "Black", NA))

gss$rur_urb_binary <- ifelse(gss$rur_urb=="Rural","Rural",
                             ifelse(gss$rur_urb=="Suburban"|gss$rur_urb=="Urban", "Non-Rural", NA))

gss$rur_urb_binary <- factor(gss$rur_urb_binary, levels = c("Rural", "Non-Rural"))

gss$religiosity <- ifelse(gss$attend=="Never"|gss$attend=="Once a year or less"|
                            gss$attend=="A few times a year"|gss$attend=="A few times a month", "Not Religious",
                          ifelse(gss$attend=="Once a week"|gss$attend=="More than once a week", "Religious", NA))

gss$religiosity <- factor(gss$religiosity, levels = c("Religious", "Not Religious"))

gss$education <- factor(gss$education, levels = c("No College Degree", "College Degree"))

# Figure 1
crosstab(gss, consci_cat, year, pct_type="col", weight=wtssall, format = "long") %>%
ggplot(aes(x = year, y = pct/100, col = consci_cat, group = consci_cat)) +  
  geom_smooth() + geom_point() +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + 
  theme_bw(base_size = 14) + xlab("") + ylab("Percentage of Americans") +
  ggtitle("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
  labs(color='Confidence in\nScientific\nCommunity') +
  scale_color_manual(values = c("gray35", "gray70", "gray0"))

# Table 1
# Calculate mean trust in all institutions by year
gss$conarmy <- car::recode(gss$conarmy, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conbus <- car::recode(gss$conbus, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conclerg <- car::recode(gss$conclerg, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$coneduc <- car::recode(gss$coneduc, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confed <- car::recode(gss$confed, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confinan <- car::recode(gss$confinan, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conjudge <- car::recode(gss$conjudge, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conlabor <- car::recode(gss$conlabor, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conlegis <- car::recode(gss$conlegis, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$contv <- car::recode(gss$contv, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confed <- car::recode(gss$confed, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conmedic <- car::recode(gss$conmedic, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conpress <- car::recode(gss$conpress, '-100=NA; -98=NA; -99=NA; -97=NA')

gss$conarmy_binary <- ifelse(gss$conarmy==1,1,0)
gss$conbus_binary <- ifelse(gss$conbus==1,1,0)
gss$conclerg_binary <- ifelse(gss$conclerg==1,1,0)
gss$coneduc_binary <- ifelse(gss$coneduc==1,1,0)
gss$confed_binary <- ifelse(gss$confed==1,1,0)
gss$confinan_binary <- ifelse(gss$confinan==1,1,0)
gss$conjudge_binary <- ifelse(gss$conjudge==1,1,0)
gss$conlabor_binary <- ifelse(gss$conlabor==1,1,0)
gss$conlegis_binary <- ifelse(gss$conlegis==1,1,0)
gss$contv_binary <- ifelse(gss$contv==1,1,0)
gss$conmedic_binary <- ifelse(gss$conmedic==1,1,0)
gss$conpress_binary <- ifelse(gss$conpress==1,1,0)

print(crosstab(gss, conarmy_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, consci_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conbus_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conclerg_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, coneduc_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, confed_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, confinan_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conjudge_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conlabor_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conlegis_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conmedic_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, conpress_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)
print(crosstab(gss, contv_binary, year, pct_type="col", weight=wtssall, format = "long"), n=68)


# Figure 2
gss$consci_gender <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Gender", 0)
gss$consci_race <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Race", 0)
gss$consci_urb <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Rurality", 0)
gss$consci_religios <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religiosity", 0)
gss$consci_educ <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Education", 0)
gss$consci_income <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Income Distribution Percentile", 0)
gss$consci_religion <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religion", 0)
gss$consci_parents_educ <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Parents' Education", 0)
gss$consci_class <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Class", 0)
gss$consci_norelig <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religion: None", 0)


p1a <- crosstab_3way(df = gss, x = sex, y = consci_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = sex, col = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2a <- crosstab_3way(df = gss, x = race_binary, y = consci_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = race_binary, col = race_binary)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3a <- crosstab_3way(df = gss, x = rur_urb_binary, y = consci_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = rur_urb_binary, col = rur_urb_binary)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4a <- crosstab_3way(df = gss, x = religiosity, y = consci_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = religiosity, col = religiosity)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5a <- crosstab_3way(df = gss, x = education, y = consci_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = education, col = education)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6a <- crosstab_3way(df = gss, x = class2, y = consci_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = class2, col = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Fig2 <- ggarrange(p1a, p2a, p3a, p4a, p5a, p6a, nrow = 2, ncol = 3, align = "v")
Fig2

# Figure 3
gss$democrat_cat <- ifelse(gss$democrat==1, "Percent Democrat", "Not Democrat")

gss$dems_gender <- ifelse(gss$democrat_cat=="Percent Democrat", "Gender", 0)
gss$dems_race <- ifelse(gss$democrat_cat=="Percent Democrat", "Race", 0)
gss$dems_urb <- ifelse(gss$democrat_cat=="Percent Democrat", "Rurality", 0)
gss$dems_religios <- ifelse(gss$democrat_cat=="Percent Democrat", "Religiosity", 0)
gss$dems_educ <- ifelse(gss$democrat_cat=="Percent Democrat", "Education", 0)
gss$dems_class <- ifelse(gss$democrat_cat=="Percent Democrat", "Class", 0)


p1b <- crosstab_3way(df = gss, x = sex, y = dems_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b <- crosstab_3way(df = gss, x = race_binary, y = dems_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b <- crosstab_3way(df = gss, x = rur_urb_binary, y = dems_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b <- crosstab_3way(df = gss, x = religiosity, y = dems_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b <- crosstab_3way(df = gss, x = education, y = dems_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(dems_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b <- crosstab_3way(df = gss, x = class2, y = dems_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Fig3 <- ggarrange(p1b, p2b, p3b, p4b, p5b, p6b, nrow = 2, ncol = 3, align = "v")
Fig3

# Figure 4
Fig4 <- crosstab_3way(df = gss, x = pid, y = consci_cat, z = year, weight = wtssall, format = "long") %>%
  ggplot(aes(year, pct/100, col = pid, linetype = pid)) +
  geom_smooth() + 
  scale_color_grey(start = 0.05, end = 0.6) +
  facet_wrap(facets = vars(consci_cat)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + ggtitle("") +
  theme(legend.position = "bottom")
Fig4

# Figure 5
#scientists
demsci <- crosstab_3way(df = gss, x = democrat, y = consci_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demsci <- demsci[, c("year", "pct")]
colnames(demsci)[2] <- "dempct"

repsci <- crosstab_3way(df = gss, x = republican, y = consci_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repsci <- repsci[, c("year", "pct")]
colnames(repsci)[2] <- "reppct"

scipol <- merge(demsci, repsci, by = "year")
scipol$gap <- scipol$dempct - scipol$reppct

p_scipol <- ggplot(data = scipol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Scientific community") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# military
demarmy <- crosstab_3way(df = gss, x = democrat, y = conarmy_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demarmy <- demarmy[, c("year", "pct")]
colnames(demarmy)[2] <- "dempct"

reparmy <- crosstab_3way(df = gss, x = republican, y = conarmy_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reparmy <- reparmy[, c("year", "pct")]
colnames(reparmy)[2] <- "reppct"

armypol <- merge(demarmy, reparmy, by = "year")
armypol$gap <- armypol$dempct - armypol$reppct

p_armypol <- ggplot(data = armypol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Military") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# major companies
dembus <- crosstab_3way(df = gss, x = democrat, y = conbus_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
dembus <- dembus[, c("year", "pct")]
colnames(dembus)[2] <- "dempct"

repbus <- crosstab_3way(df = gss, x = republican, y = conbus_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repbus <- repbus[, c("year", "pct")]
colnames(repbus)[2] <- "reppct"

buspol <- merge(dembus, repbus, by = "year")
buspol$gap <- buspol$dempct - buspol$reppct

p_buspol <- ggplot(data = buspol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Major companies") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized religion
demclerg <- crosstab_3way(df = gss, x = democrat, y = conclerg_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demclerg <- demclerg[, c("year", "pct")]
colnames(demclerg)[2] <- "dempct"

repclerg <- crosstab_3way(df = gss, x = republican, y = conclerg_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repclerg <- repclerg[, c("year", "pct")]
colnames(repclerg)[2] <- "reppct"

clergpol <- merge(demclerg, repclerg, by = "year")
clergpol$gap <- clergpol$dempct - clergpol$reppct

p_clergpol <- ggplot(data = clergpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized religion") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# education
demeduc <- crosstab_3way(df = gss, x = democrat, y = coneduc_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demeduc <- demeduc[, c("year", "pct")]
colnames(demeduc)[2] <- "dempct"

repeduc <- crosstab_3way(df = gss, x = republican, y = coneduc_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repeduc <- repeduc[, c("year", "pct")]
colnames(repeduc)[2] <- "reppct"

educpol <- merge(demeduc, repeduc, by = "year")
educpol$gap <- educpol$dempct - educpol$reppct

p_educpol <- ggplot(data = educpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Education") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# banks and financial institutions
demfinan <- crosstab_3way(df = gss, x = democrat, y = confinan_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demfinan <- demfinan[, c("year", "pct")]
colnames(demfinan)[2] <- "dempct"

repfinan <- crosstab_3way(df = gss, x = republican, y = confinan_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repfinan <- repfinan[, c("year", "pct")]
colnames(repfinan)[2] <- "reppct"

finanpol <- merge(demfinan, repfinan, by = "year")
finanpol$gap <- finanpol$dempct - finanpol$reppct

p_finanpol <- ggplot(data = finanpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Banks and financial institutions") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# supreme court
demjudge <- crosstab_3way(df = gss, x = democrat, y = conjudge_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demjudge <- demjudge[, c("year", "pct")]
colnames(demjudge)[2] <- "dempct"

repjudge <- crosstab_3way(df = gss, x = republican, y = conjudge_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repjudge <- repjudge[, c("year", "pct")]
colnames(repjudge)[2] <- "reppct"

judgepol <- merge(demjudge, repjudge, by = "year")
judgepol$gap <- judgepol$dempct - judgepol$reppct

p_judgepol <- ggplot(data = judgepol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Supreme Court") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized labor
demlabor <- crosstab_3way(df = gss, x = democrat, y = conlabor_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demlabor <- demlabor[, c("year", "pct")]
colnames(demlabor)[2] <- "dempct"

replabor <- crosstab_3way(df = gss, x = republican, y = conlabor_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
replabor <- replabor[, c("year", "pct")]
colnames(replabor)[2] <- "reppct"

laborpol <- merge(demlabor, replabor, by = "year")
laborpol$gap <- laborpol$dempct - laborpol$reppct

p_laborpol <- ggplot(data = laborpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized labor") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Congress
demlegis <- crosstab_3way(df = gss, x = democrat, y = conlegis_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demlegis <- demlegis[, c("year", "pct")]
colnames(demlegis)[2] <- "dempct"

replegis <- crosstab_3way(df = gss, x = republican, y = conlegis_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
replegis <- replegis[, c("year", "pct")]
colnames(replegis)[2] <- "reppct"

legispol <- merge(demlegis, replegis, by = "year")
legispol$gap <- legispol$dempct - legispol$reppct

p_legispol <- ggplot(data = legispol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Congress") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Medicine
demmedic <- crosstab_3way(df = gss, x = democrat, y = conmedic_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demmedic <- demmedic[, c("year", "pct")]
colnames(demmedic)[2] <- "dempct"

repmedic <- crosstab_3way(df = gss, x = republican, y = conmedic_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repmedic <- repmedic[, c("year", "pct")]
colnames(repmedic)[2] <- "reppct"

medicpol <- merge(demmedic, repmedic, by = "year")
medicpol$gap <- medicpol$dempct - medicpol$reppct

p_medicpol <- ggplot(data = medicpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Medicine") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# press
dempress <- crosstab_3way(df = gss, x = democrat, y = conpress_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
dempress <- dempress[, c("year", "pct")]
colnames(dempress)[2] <- "dempct"

reppress <- crosstab_3way(df = gss, x = republican, y = conpress_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reppress <- reppress[, c("year", "pct")]
colnames(reppress)[2] <- "reppct"

presspol <- merge(dempress, reppress, by = "year")
presspol$gap <- presspol$dempct - presspol$reppct

p_presspol <- ggplot(data = presspol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Press") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# television
demtv <- crosstab_3way(df = gss, x = democrat, y = contv_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demtv <- demtv[, c("year", "pct")]
colnames(demtv)[2] <- "dempct"

reptv <- crosstab_3way(df = gss, x = republican, y = contv_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reptv <- reptv[, c("year", "pct")]
colnames(reptv)[2] <- "reppct"

tvpol <- merge(demtv, reptv, by = "year")
tvpol$gap <- tvpol$dempct - tvpol$reppct

p_tvpol <- ggplot(data = tvpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Television") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

Fig5 <- ggarrange(p_scipol, p_educpol, p_presspol, p_tvpol,
                      p_judgepol, p_medicpol, p_armypol, p_clergpol,
                      p_buspol, p_laborpol, p_finanpol, p_legispol,
                      nrow = 3, ncol = 4, align = "v")
Fig5 <- annotate_figure(Fig5,
                            left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))
Fig5

# Figure 6
fig6data <- read_dta("fig6data.dta")

scitrustmeans <- c(mean(fig6data$confsci[fig6data$rep==1]), mean(fig6data$confsci[fig6data$dem==1]),
                   mean(fig6data$pconfsci[fig6data$rep==1]), mean(fig6data$pconfsci[fig6data$dem==1]),
                   mean(fig6data$scidiff[fig6data$rep==1]), mean(fig6data$scidiff[fig6data$dem==1]))
parties <- c("Republicans", "Democrats", "Republicans", "Democrats", "Republicans", "Democrats")
cat <- c("Trust in Scientists", "Trust in Scientists", "Perception of Out-Party Trust in Scientists", "Perception of Out-Party Trust in Scientists", "Trust Differential", "Trust Differential")
sds <- c(sd(fig6data$confsci[fig6data$rep==1]), sd(fig6data$confsci[fig6data$dem==1]),
         sd(fig6data$pconfsci[fig6data$rep==1]), sd(fig6data$pconfsci[fig6data$dem==1]),
         sd(fig6data$scidiff[fig6data$rep==1]), sd(fig6data$scidiff[fig6data$dem==1]))
n <- c(742, 759, 742, 759, 742, 759)

df <- data.frame(scitrustmeans,parties,cat, sds, n)
df$se <- qnorm(0.975) * (df$sds/sqrt(df$n))
df$Lower_CI <- (df$scitrustmeans - df$se)
df$Upper_CI <- (df$scitrustmeans + df$se)

df$cat <- factor(df$cat, levels = c("Trust Differential", "Perception of Out-Party Trust in Scientists", "Trust in Scientists"))

Fig6 <- ggplot(df, aes(parties, scitrustmeans)) +
  geom_col(aes(fill = forcats::fct_rev(cat)), position = position_dodge(0.7), 
           color = "black", alpha = 0.7, width = 0.6) +
  geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI, group = forcats::fct_rev(cat)), 
                position = position_dodge(0.7), width = 0.2) +
  scale_fill_manual(values = c("black", "gray30", "gray70")) +
  theme_minimal(base_size = 20) +
  theme(panel.grid.major.x = element_blank()) +
  labs(x = NULL, y = "Average Score (0 to 100 scale)") + ylim(0,100) +
  theme(legend.position="bottom") +
  theme(legend.title=element_blank()) +
  annotate("text", x = 0.765, y = 86.34, label = "76.34", color="black", size=6) +
  annotate("text", x = 1, y = 45.84, label = "35.84", color="black", size=6) +
  annotate("text", x = 1.235, y = 53.34, label = "43.34", color="black", size=6) +
  annotate("text", x = 1.765, y = 65.29, label = "55.29", color="black", size=6) +
  annotate("text", x = 2, y = 65.74, label = "55.74", color="black", size=6) +
  annotate("text", x = 2.235, y = 36.64, label = "26.64", color="black", size=6)

Fig6


# Figure 7. Regression Coefficients with Trust in Scientists and Researchers.
fig7data <- read_stata("fig7data")
fig7data_reps <- subset(fig7data, fig7data$rep==1)
fig7data_dems <- subset(fig7data, fig7data$dem==1)

fig7data$pid=case_when(
  fig7data$rep==1~"Republican",
  fig7data$dem==1~"Democrat",
)

fig7_ovr <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data)
fig7_dems <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data_dems)
fig7_reps <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data_reps)

Fig7 <- plot_summs(fig7_ovr,fig7_dems,fig7_reps,
           model.names=c("Overall","Democrats","Republicans"),
           colors = c("black", "gray40", "gray70"),
           coefs = c("Medical Scientists"="medicaltr",
                     "Applied Scientists"="applied",
                     "Natural Scientists"="basic",
                     "Social Scientists"="socialsci",
                     "Climate Scientists"="climatetr",
                     "COVID-19 Biochemists"="covidtr",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.2,0.6) + xlab("Coefficients") + ylab("Trust Target") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 16), axis.text.y = element_text(size = 16))

Fig7


# Table 2
tab2data <- read_dta("tab2data")

tab2data$naturalim <- rowMeans(subset(tab2data, select=c(physicistim, chemim, biologim)), na.rm = TRUE)

mean(tab2data$scimethodim)
sd(tab2data$scimethodim)
mean(tab2data$scimethodim[tab2data$pidrep<3.5])
sd(tab2data$scimethodim[tab2data$pidrep<3.5])
mean(tab2data$scimethodim[tab2data$pidrep>4.5])
sd(tab2data$scimethodim[tab2data$pidrep>4.5])


mean(tab2data$doctorim)
sd(tab2data$doctorim)
mean(tab2data$doctorim[tab2data$pidrep<3.5])
sd(tab2data$doctorim[tab2data$pidrep<3.5])
mean(tab2data$doctorim[tab2data$pidrep>4.5])
sd(tab2data$doctorim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$doctorim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$doctorim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$doctorim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$engineerim)
sd(tab2data$engineerim)
mean(tab2data$engineerim[tab2data$pidrep<3.5])
sd(tab2data$engineerim[tab2data$pidrep<3.5])
mean(tab2data$engineerim[tab2data$pidrep>4.5])
sd(tab2data$engineerim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$engineerim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$engineerim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$engineerim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$naturalim)
sd(tab2data$naturalim)
mean(tab2data$naturalim[tab2data$pidrep<3.5])
sd(tab2data$naturalim[tab2data$pidrep<3.5])
mean(tab2data$naturalim[tab2data$pidrep>4.5])
sd(tab2data$naturalim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$naturalim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$naturalim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$naturalim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$climateim)
sd(tab2data$climateim)
mean(tab2data$climateim[tab2data$pidrep<3.5])
sd(tab2data$climateim[tab2data$pidrep<3.5])
mean(tab2data$climateim[tab2data$pidrep>4.5])
sd(tab2data$climateim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$climateim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$climateim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$climateim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$covidim)
sd(tab2data$covidim)
mean(tab2data$covidim[tab2data$pidrep<3.5])
sd(tab2data$covidim[tab2data$pidrep<3.5])
mean(tab2data$covidim[tab2data$pidrep>4.5])
sd(tab2data$covidim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$covidim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$covidim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$covidim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$professorim)
sd(tab2data$professorim)
mean(tab2data$professorim[tab2data$pidrep<3.5])
sd(tab2data$professorim[tab2data$pidrep<3.5])
mean(tab2data$professorim[tab2data$pidrep>4.5])
sd(tab2data$professorim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$professorim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$professorim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$professorim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$highschoolim)
sd(tab2data$highschoolim)
mean(tab2data$highschoolim[tab2data$pidrep<3.5])
sd(tab2data$highschoolim[tab2data$pidrep<3.5])
mean(tab2data$highschoolim[tab2data$pidrep>4.5])
sd(tab2data$highschoolim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$highschoolim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$highschoolim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$highschoolim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$popularim)
sd(tab2data$popularim)
mean(tab2data$popularim[tab2data$pidrep<3.5])
sd(tab2data$popularim[tab2data$pidrep<3.5])
mean(tab2data$popularim[tab2data$pidrep>4.5])
sd(tab2data$popularim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$popularim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$popularim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$popularim[tab2data$pidrep>4.5], alternative = "two.sided")

################
### Appendix ###
################

# Looking for Mode Effects
print(crosstab_3way(df = gss, x = mode, y = consci_binary, z = year, weight = wtssall, format = "long", remove = 0), n=98)

# Figure 1
# Predictors of confidence in scientists by decade
gss1970s <- subset(gss, gss$year<1979.5)
gss1980s <- subset(gss, gss$year>1979.5 & gss$year<1989.5)
gss1990s <- subset(gss, gss$year>1989.5 & gss$year<1999.5)
gss2000s <- subset(gss, gss$year>1999.5 & gss$year<2009.5)
gss2010s <- subset(gss, gss$year>2009.5 & gss$year<2019.5)

sci_70s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
sci_80s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
sci_90s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
sci_00s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
sci_10s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(sci_70s,sci_80s,sci_90s,sci_00s,sci_10s,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.1, 0.25) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))

# Figure 2
# Looking at lowest confidence in science (7-20-24)
gss$consci_gender2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Gender", 0)
gss$consci_race2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Race", 0)
gss$consci_urb2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Rurality", 0)
gss$consci_religios2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religiosity", 0)
gss$consci_educ2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Education", 0)
gss$consci_income2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Income Distribution Percentile", 0)
gss$consci_religion2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religion", 0)
gss$consci_parents_educ2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Parents' Education", 0)
gss$consci_class2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Class", 0)
gss$consci_norelig2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religion: None", 0)


p1b <- crosstab_3way(df = gss, x = sex, y = consci_gender2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_gender2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b <- crosstab_3way(df = gss, x = race_binary, y = consci_race2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_race2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b <- crosstab_3way(df = gss, x = rur_urb_binary, y = consci_urb2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_urb2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b <- crosstab_3way(df = gss, x = religiosity, y = consci_religios2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_religios2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b <- crosstab_3way(df = gss, x = education, y = consci_educ2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(consci_educ2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b <- crosstab_3way(df = gss, x = class2, y = consci_class2, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_class2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_dem_lowest <- ggarrange(p1b, p2b, p3b, p4b, p5b, p6b, nrow = 2, ncol = 3, align = "v")
fig_dem_lowest


# Percent Identifying as Republicans by Party over Time
gss$republican_cat <- ifelse(gss$republican==1, "Percent republican", "Not republican")

gss$reps_gender <- ifelse(gss$republican_cat=="Percent republican", "Gender", 0)
gss$reps_race <- ifelse(gss$republican_cat=="Percent republican", "Race", 0)
gss$reps_urb <- ifelse(gss$republican_cat=="Percent republican", "Rurality", 0)
gss$reps_religios <- ifelse(gss$republican_cat=="Percent republican", "Religiosity", 0)
gss$reps_educ <- ifelse(gss$republican_cat=="Percent republican", "Education", 0)
gss$reps_class <- ifelse(gss$republican_cat=="Percent republican", "Class", 0)


p1c <- crosstab_3way(df = gss, x = sex, y = reps_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2c <- crosstab_3way(df = gss, x = race_binary, y = reps_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3c <- crosstab_3way(df = gss, x = rur_urb_binary, y = reps_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4c <- crosstab_3way(df = gss, x = religiosity, y = reps_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5c <- crosstab_3way(df = gss, x = education, y = reps_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(reps_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6c <- crosstab_3way(df = gss, x = class2, y = reps_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_republican <- ggarrange(p1c, p2c, p3c, p4c, p5c, p6c, nrow = 2, ncol = 3, align = "v")
fig_republican <- annotate_figure(fig_republican,
                                  top=text_grob("Percent Republican by Group", size = 18))
fig_republican


# Percent Identifying as Independent over Time
gss$independent_cat <- ifelse(gss$independent==1, "Percent independent", "Not independent")

gss$inds_gender <- ifelse(gss$independent_cat=="Percent independent", "Gender", 0)
gss$inds_race <- ifelse(gss$independent_cat=="Percent independent", "Race", 0)
gss$inds_urb <- ifelse(gss$independent_cat=="Percent independent", "Rurality", 0)
gss$inds_religios <- ifelse(gss$independent_cat=="Percent independent", "Religiosity", 0)
gss$inds_educ <- ifelse(gss$independent_cat=="Percent independent", "Education", 0)
gss$inds_class <- ifelse(gss$independent_cat=="Percent independent", "Class", 0)


p1d <- crosstab_3way(df = gss, x = sex, y = inds_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2d <- crosstab_3way(df = gss, x = race_binary, y = inds_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3d <- crosstab_3way(df = gss, x = rur_urb_binary, y = inds_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4d <- crosstab_3way(df = gss, x = religiosity, y = inds_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5d <- crosstab_3way(df = gss, x = education, y = inds_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(inds_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6d <- crosstab_3way(df = gss, x = class2, y = inds_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_independent <- ggarrange(p1d, p2d, p3d, p4d, p5d, p6d, nrow = 2, ncol = 3, align = "v")
fig_independent <- annotate_figure(fig_independent,
                                  top=text_grob("Percent Independent by Group", size = 18))
fig_independent

# Predictors of being a Democrat by decade
dem_70s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
dem_80s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
dem_90s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
dem_00s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
dem_10s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(dem_70s,dem_80s,dem_90s,dem_00s,dem_10s,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.4, 0.2) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))


# Results With Independent Leaners Treated as Independents Rather than Partisans 
# Democrats
gss$democrat_noleaners_cat <- ifelse(gss$democrat_noleaners==1, "Percent Democrat", "Not Democrat")

gss$dems_noleaners_gender <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Gender", 0)
gss$dems_noleaners_race <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Race", 0)
gss$dems_noleaners_urb <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Rurality", 0)
gss$dems_noleaners_religios <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Religiosity", 0)
gss$dems_noleaners_educ <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Education", 0)
gss$dems_noleaners_class <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Class", 0)


p1a_app <- crosstab_3way(df = gss, x = sex, y = dems_noleaners_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2a_app <- crosstab_3way(df = gss, x = race_binary, y = dems_noleaners_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3a_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = dems_noleaners_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4a_app <- crosstab_3way(df = gss, x = religiosity, y = dems_noleaners_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5a_app <- crosstab_3way(df = gss, x = education, y = dems_noleaners_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(dems_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6a_app <- crosstab_3way(df = gss, x = class2, y = dems_noleaners_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1a_app, p2a_app, p3a_app, p4a_app, p5a_app, p6a_app, nrow = 2, ncol = 3, align = "v")


# Republicans
gss$republican_noleaners_cat <- ifelse(gss$republican_noleaners==1, "Percent Republican", "Not Republican")

gss$reps_noleaners_gender <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Gender", 0)
gss$reps_noleaners_race <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Race", 0)
gss$reps_noleaners_urb <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Rurality", 0)
gss$reps_noleaners_religios <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Religiosity", 0)
gss$reps_noleaners_educ <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Education", 0)
gss$reps_noleaners_class <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Class", 0)


p1b_app <- crosstab_3way(df = gss, x = sex, y = reps_noleaners_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b_app <- crosstab_3way(df = gss, x = race_binary, y = reps_noleaners_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = reps_noleaners_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b_app <- crosstab_3way(df = gss, x = religiosity, y = reps_noleaners_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b_app <- crosstab_3way(df = gss, x = education, y = reps_noleaners_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(reps_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b_app <- crosstab_3way(df = gss, x = class2, y = reps_noleaners_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1b_app, p2b_app, p3b_app, p4b_app, p5b_app, p6b_app, nrow = 2, ncol = 3, align = "v")


# independents
gss$independent_noleaners_cat <- ifelse(gss$independent_noleaners==1, "Percent Independent", "Not independent")

gss$inds_noleaners_gender <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Gender", 0)
gss$inds_noleaners_race <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Race", 0)
gss$inds_noleaners_urb <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Rurality", 0)
gss$inds_noleaners_religios <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Religiosity", 0)
gss$inds_noleaners_educ <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Education", 0)
gss$inds_noleaners_class <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Class", 0)


p1c_app <- crosstab_3way(df = gss, x = sex, y = inds_noleaners_gender, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2c_app <- crosstab_3way(df = gss, x = race_binary, y = inds_noleaners_race, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3c_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = inds_noleaners_urb, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4c_app <- crosstab_3way(df = gss, x = religiosity, y = inds_noleaners_religios, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5c_app <- crosstab_3way(df = gss, x = education, y = inds_noleaners_educ, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(inds_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6c_app <- crosstab_3way(df = gss, x = class2, y = inds_noleaners_class, z = year, weight = wtssall, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1c_app, p2c_app, p3c_app, p4c_app, p5c_app, p6c_app, nrow = 2, ncol = 3, align = "v")

# Regression of Democrat Affiliation (Without Leaners as Partisans) By Decade 
gss1970s <- subset(gss, gss$year<1979.5)
gss1980s <- subset(gss, gss$year>1979.5 & gss$year<1989.5)
gss1990s <- subset(gss, gss$year>1989.5 & gss$year<1999.5)
gss2000s <- subset(gss, gss$year>1999.5 & gss$year<2009.5)
gss2010s <- subset(gss, gss$year>2009.5 & gss$year<2018.5)

dem_70s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
dem_80s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
dem_90s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
dem_00s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
dem_10s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(dem_70s_app,dem_80s_app,dem_90s_app,dem_00s_app,dem_10s_app,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.4, 0.2) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))


# Confidence in the Scientific Community By Party (Without Leaners as Partisans) Over Time
crosstab_3way(df = gss, x = partyid_without_leaners_as_partisans, y = consci_cat, z = year, weight = wtssall, format = "long") %>%
  ggplot(aes(year, pct/100, col = partyid_without_leaners_as_partisans, linetype = partyid_without_leaners_as_partisans)) + 
  geom_smooth() + 
  scale_color_grey(start = 0.05, end = 0.6) +
  facet_wrap(facets = vars(consci_cat)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + ggtitle("") +
  theme(legend.position = "bottom")


# Confidence in Scientists by Party (without leaners as partisans)
df <- crosstab_3way(df = gss, x = partyid_without_leaners_as_partisans, y = consci_binary, z = year, weight = wtssall, remove=0, format = "long")

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2018])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2018])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2021])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2021])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2022])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2022])

# Partisan Gaps (Without Leaners as Partisans) in Confidence in Institutions over Time
#scientists
demsci <- crosstab_3way(df = gss, x = democrat_noleaners, y = consci_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demsci <- demsci[, c("year", "pct")]
colnames(demsci)[2] <- "dempct"

repsci <- crosstab_3way(df = gss, x = republican_noleaners, y = consci_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repsci <- repsci[, c("year", "pct")]
colnames(repsci)[2] <- "reppct"

scipol <- merge(demsci, repsci, by = "year")
scipol$gap <- scipol$dempct - scipol$reppct

p_scipol <- ggplot(data = scipol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Scientific community") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# military
demarmy <- crosstab_3way(df = gss, x = democrat_noleaners, y = conarmy_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demarmy <- demarmy[, c("year", "pct")]
colnames(demarmy)[2] <- "dempct"

reparmy <- crosstab_3way(df = gss, x = republican_noleaners, y = conarmy_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reparmy <- reparmy[, c("year", "pct")]
colnames(reparmy)[2] <- "reppct"

armypol <- merge(demarmy, reparmy, by = "year")
armypol$gap <- armypol$dempct - armypol$reppct

p_armypol <- ggplot(data = armypol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Military") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# major companies
dembus <- crosstab_3way(df = gss, x = democrat_noleaners, y = conbus_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
dembus <- dembus[, c("year", "pct")]
colnames(dembus)[2] <- "dempct"

repbus <- crosstab_3way(df = gss, x = republican_noleaners, y = conbus_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repbus <- repbus[, c("year", "pct")]
colnames(repbus)[2] <- "reppct"

buspol <- merge(dembus, repbus, by = "year")
buspol$gap <- buspol$dempct - buspol$reppct

p_buspol <- ggplot(data = buspol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Major companies") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized religion
demclerg <- crosstab_3way(df = gss, x = democrat_noleaners, y = conclerg_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demclerg <- demclerg[, c("year", "pct")]
colnames(demclerg)[2] <- "dempct"

repclerg <- crosstab_3way(df = gss, x = republican_noleaners, y = conclerg_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repclerg <- repclerg[, c("year", "pct")]
colnames(repclerg)[2] <- "reppct"

clergpol <- merge(demclerg, repclerg, by = "year")
clergpol$gap <- clergpol$dempct - clergpol$reppct

p_clergpol <- ggplot(data = clergpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized religion") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# education
demeduc <- crosstab_3way(df = gss, x = democrat_noleaners, y = coneduc_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demeduc <- demeduc[, c("year", "pct")]
colnames(demeduc)[2] <- "dempct"

repeduc <- crosstab_3way(df = gss, x = republican_noleaners, y = coneduc_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repeduc <- repeduc[, c("year", "pct")]
colnames(repeduc)[2] <- "reppct"

educpol <- merge(demeduc, repeduc, by = "year")
educpol$gap <- educpol$dempct - educpol$reppct

p_educpol <- ggplot(data = educpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Education") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# executive branch
demfed <- crosstab_3way(df = gss, x = democrat_noleaners, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican_noleaners, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Partisan Gap in Executive Branch") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

p_fedpol <- annotate_figure(p_fedpol,
                            left = text_grob("democrat_noleanerss with Great Deal of Confidence -\nrepublican_noleanerss with Great Deal of Confidence", size = 14, rot = 90))

# banks and financial institutions
demfinan <- crosstab_3way(df = gss, x = democrat_noleaners, y = confinan_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demfinan <- demfinan[, c("year", "pct")]
colnames(demfinan)[2] <- "dempct"

repfinan <- crosstab_3way(df = gss, x = republican_noleaners, y = confinan_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repfinan <- repfinan[, c("year", "pct")]
colnames(repfinan)[2] <- "reppct"

finanpol <- merge(demfinan, repfinan, by = "year")
finanpol$gap <- finanpol$dempct - finanpol$reppct

p_finanpol <- ggplot(data = finanpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Banks and financial institutions") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# supreme court
demjudge <- crosstab_3way(df = gss, x = democrat_noleaners, y = conjudge_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demjudge <- demjudge[, c("year", "pct")]
colnames(demjudge)[2] <- "dempct"

repjudge <- crosstab_3way(df = gss, x = republican_noleaners, y = conjudge_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repjudge <- repjudge[, c("year", "pct")]
colnames(repjudge)[2] <- "reppct"

judgepol <- merge(demjudge, repjudge, by = "year")
judgepol$gap <- judgepol$dempct - judgepol$reppct

p_judgepol <- ggplot(data = judgepol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Supreme Court") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized labor
demlabor <- crosstab_3way(df = gss, x = democrat_noleaners, y = conlabor_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demlabor <- demlabor[, c("year", "pct")]
colnames(demlabor)[2] <- "dempct"

replabor <- crosstab_3way(df = gss, x = republican_noleaners, y = conlabor_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
replabor <- replabor[, c("year", "pct")]
colnames(replabor)[2] <- "reppct"

laborpol <- merge(demlabor, replabor, by = "year")
laborpol$gap <- laborpol$dempct - laborpol$reppct

p_laborpol <- ggplot(data = laborpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized labor") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Congress
demlegis <- crosstab_3way(df = gss, x = democrat_noleaners, y = conlegis_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demlegis <- demlegis[, c("year", "pct")]
colnames(demlegis)[2] <- "dempct"

replegis <- crosstab_3way(df = gss, x = republican_noleaners, y = conlegis_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
replegis <- replegis[, c("year", "pct")]
colnames(replegis)[2] <- "reppct"

legispol <- merge(demlegis, replegis, by = "year")
legispol$gap <- legispol$dempct - legispol$reppct

p_legispol <- ggplot(data = legispol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Congress") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Medicine
demmedic <- crosstab_3way(df = gss, x = democrat_noleaners, y = conmedic_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demmedic <- demmedic[, c("year", "pct")]
colnames(demmedic)[2] <- "dempct"

repmedic <- crosstab_3way(df = gss, x = republican_noleaners, y = conmedic_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repmedic <- repmedic[, c("year", "pct")]
colnames(repmedic)[2] <- "reppct"

medicpol <- merge(demmedic, repmedic, by = "year")
medicpol$gap <- medicpol$dempct - medicpol$reppct

p_medicpol <- ggplot(data = medicpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Medicine") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# press
dempress <- crosstab_3way(df = gss, x = democrat_noleaners, y = conpress_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
dempress <- dempress[, c("year", "pct")]
colnames(dempress)[2] <- "dempct"

reppress <- crosstab_3way(df = gss, x = republican_noleaners, y = conpress_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reppress <- reppress[, c("year", "pct")]
colnames(reppress)[2] <- "reppct"

presspol <- merge(dempress, reppress, by = "year")
presspol$gap <- presspol$dempct - presspol$reppct

p_presspol <- ggplot(data = presspol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Press") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# television
demtv <- crosstab_3way(df = gss, x = democrat_noleaners, y = contv_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demtv <- demtv[, c("year", "pct")]
colnames(demtv)[2] <- "dempct"

reptv <- crosstab_3way(df = gss, x = republican_noleaners, y = contv_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
reptv <- reptv[, c("year", "pct")]
colnames(reptv)[2] <- "reppct"

tvpol <- merge(demtv, reptv, by = "year")
tvpol$gap <- tvpol$dempct - tvpol$reppct

p_tvpol <- ggplot(data = tvpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Television") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

FigPartisanGap_app <- ggarrange(p_scipol, p_educpol, p_presspol, p_tvpol,
                  p_judgepol, p_medicpol, p_armypol, p_clergpol,
                  p_buspol, p_laborpol, p_finanpol, p_legispol,
                  nrow = 3, ncol = 4, align = "v")
FigPartisanGap_app <- annotate_figure(FigPartisanGap_app,
                        left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))
FigPartisanGap_app

# Executive branch plots
demfed <- crosstab_3way(df = gss, x = democrat, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

annotate_figure(p_fedpol,
                left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))

# Executive branch plot without leaners as partisans 
demfed <- crosstab_3way(df = gss, x = democrat_noleaners, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican_noleaners, y = confed_binary, z = year, weight = wtssall, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

annotate_figure(p_fedpol,
                left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))



# Perceptions of Institutions Data

# Trust Differentials table in appendix
# scientific community
mean(fig6data$scidiff)
sd(fig6data$scidiff)
mean(fig6data$scidiff[fig6data$dem==1])
sd(fig6data$scidiff[fig6data$dem==1])
mean(fig6data$scidiff[fig6data$rep==1])
sd(fig6data$scidiff[fig6data$rep==1])

# 
mean(fig6data$educdiff)
sd(fig6data$educdiff)
mean(fig6data$educdiff[fig6data$dem==1])
sd(fig6data$educdiff[fig6data$dem==1])
mean(fig6data$educdiff[fig6data$rep==1])
sd(fig6data$educdiff[fig6data$rep==1])

# 
mean(fig6data$mediadiff)
sd(fig6data$mediadiff)
mean(fig6data$mediadiff[fig6data$dem==1])
sd(fig6data$mediadiff[fig6data$dem==1])
mean(fig6data$mediadiff[fig6data$rep==1])
sd(fig6data$mediadiff[fig6data$rep==1])

# 
mean(fig6data$socmeddiff)
sd(fig6data$socmeddiff)
mean(fig6data$socmeddiff[fig6data$dem==1])
sd(fig6data$socmeddiff[fig6data$dem==1])
mean(fig6data$socmeddiff[fig6data$rep==1])
sd(fig6data$socmeddiff[fig6data$rep==1])

# 
mean(fig6data$militarydiff)
sd(fig6data$militarydiff)
mean(fig6data$militarydiff[fig6data$dem==1])
sd(fig6data$militarydiff[fig6data$dem==1])
mean(fig6data$militarydiff[fig6data$rep==1])
sd(fig6data$militarydiff[fig6data$rep==1])

# 
mean(fig6data$policediff)
sd(fig6data$policediff)
mean(fig6data$policediff[fig6data$dem==1])
sd(fig6data$policediff[fig6data$dem==1])
mean(fig6data$policediff[fig6data$rep==1])
sd(fig6data$policediff[fig6data$rep==1])

# 
mean(fig6data$docotorsdiff)
sd(fig6data$docotorsdiff)
mean(fig6data$docotorsdiff[fig6data$dem==1])
sd(fig6data$docotorsdiff[fig6data$dem==1])
mean(fig6data$docotorsdiff[fig6data$rep==1])
sd(fig6data$docotorsdiff[fig6data$rep==1])

# 
mean(fig6data$religdiff)
sd(fig6data$religdiff)
mean(fig6data$religdiff[fig6data$dem==1])
sd(fig6data$religdiff[fig6data$dem==1])
mean(fig6data$religdiff[fig6data$rep==1])
sd(fig6data$religdiff[fig6data$rep==1])

# 
mean(fig6data$labordiff)
sd(fig6data$labordiff)
mean(fig6data$labordiff[fig6data$dem==1])
sd(fig6data$labordiff[fig6data$dem==1])
mean(fig6data$labordiff[fig6data$rep==1])
sd(fig6data$labordiff[fig6data$rep==1])

# 
mean(fig6data$companiesdiff)
sd(fig6data$companiesdiff)
mean(fig6data$companiesdiff[fig6data$dem==1])
sd(fig6data$companiesdiff[fig6data$dem==1])
mean(fig6data$companiesdiff[fig6data$rep==1])
sd(fig6data$companiesdiff[fig6data$rep==1])

# 
mean(fig6data$courtdiff)
sd(fig6data$courtdiff)
mean(fig6data$courtdiff[fig6data$dem==1])
sd(fig6data$courtdiff[fig6data$dem==1])
mean(fig6data$courtdiff[fig6data$rep==1])
sd(fig6data$courtdiff[fig6data$rep==1])

# 
mean(fig6data$congressdiff)
sd(fig6data$congressdiff)
mean(fig6data$congressdiff[fig6data$dem==1])
sd(fig6data$congressdiff[fig6data$dem==1])
mean(fig6data$congressdiff[fig6data$rep==1])
sd(fig6data$congressdiff[fig6data$rep==1])

# 
mean(fig6data$presidentdiff)
sd(fig6data$presidentdiff)
mean(fig6data$presidentdiff[fig6data$dem==1])
sd(fig6data$presidentdiff[fig6data$dem==1])
mean(fig6data$presidentdiff[fig6data$rep==1])
sd(fig6data$presidentdiff[fig6data$rep==1])

# Who are Scientists Data

# Tables Underlying Coefficient Plots
# Figure 7
summary(fig7_ovr)
summary(fig7_dems)
summary(fig7_reps)

stargazer::stargazer(fig7_ovr, fig7_dems, fig7_reps,
                     column.labels=c("All Respondents", "Democrats", "Republicans"),
                     dep.var.labels = "General Trust in Scientists",
                     covariate.labels=c("Constant", "Medical Scientists", "Applied Scientists", "Natural Scientists",
                                        "Social Scientists", "Climate Scientists", "COVID-19 Biochemists"),
                     summary = TRUE, align=TRUE,
                     type="html",
                     font.size="scriptsize", 
                     table.placement="H",
                     column.sep.width = "1pt", 
                     dep.var.labels.include = T,
                     intercept.bottom = FALSE,
                     out = "C:/Users/jonsc/Downloads/test.html", report=('vcsp'),
                     digits = 3,  no.space=TRUE)


# Trust Scores
mean(tab2data$scienctr)
sd(tab2data$scienctr)
mean(tab2data$scienctr[tab2data$pidrep<3.5])
sd(tab2data$scienctr[tab2data$pidrep<3.5])
mean(tab2data$scienctr[tab2data$pidrep>4.5])
sd(tab2data$scienctr[tab2data$pidrep>4.5])

mean(tab2data$medicaltr)
sd(tab2data$medicaltr)
mean(tab2data$medicaltr[tab2data$pidrep<3.5])
sd(tab2data$medicaltr[tab2data$pidrep<3.5])
mean(tab2data$medicaltr[tab2data$pidrep>4.5])
sd(tab2data$medicaltr[tab2data$pidrep>4.5])

t.test(tab2data$medicaltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$medicaltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$medicaltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$appliedtr <- rowMeans(subset(tab2data, select=c(aerotr, civiltr, computertr)), na.rm = TRUE)

mean(tab2data$appliedtr)
sd(tab2data$appliedtr)
mean(tab2data$appliedtr[tab2data$pidrep<3.5])
sd(tab2data$appliedtr[tab2data$pidrep<3.5])
mean(tab2data$appliedtr[tab2data$pidrep>4.5])
sd(tab2data$appliedtr[tab2data$pidrep>4.5])

t.test(tab2data$appliedtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$appliedtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$appliedtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$aerotr)
sd(tab2data$aerotr)
mean(tab2data$aerotr[tab2data$pidrep<3.5])
sd(tab2data$aerotr[tab2data$pidrep<3.5])
mean(tab2data$aerotr[tab2data$pidrep>4.5])
sd(tab2data$aerotr[tab2data$pidrep>4.5])

t.test(tab2data$aerotr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$aerotr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$aerotr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$civiltr)
sd(tab2data$civiltr)
mean(tab2data$civiltr[tab2data$pidrep<3.5])
sd(tab2data$civiltr[tab2data$pidrep<3.5])
mean(tab2data$civiltr[tab2data$pidrep>4.5])
sd(tab2data$civiltr[tab2data$pidrep>4.5])

t.test(tab2data$civiltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$civiltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$civiltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$computertr)
sd(tab2data$computertr)
mean(tab2data$computertr[tab2data$pidrep<3.5])
sd(tab2data$computertr[tab2data$pidrep<3.5])
mean(tab2data$computertr[tab2data$pidrep>4.5])
sd(tab2data$computertr[tab2data$pidrep>4.5])

t.test(tab2data$computertr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$computertr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$computertr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$naturaltr <- rowMeans(subset(tab2data, select=c(physictr, chemisttr, biologisttr)), na.rm = TRUE)

mean(tab2data$naturaltr)
sd(tab2data$naturaltr)
mean(tab2data$naturaltr[tab2data$pidrep<3.5])
sd(tab2data$naturaltr[tab2data$pidrep<3.5])
mean(tab2data$naturaltr[tab2data$pidrep>4.5])
sd(tab2data$naturaltr[tab2data$pidrep>4.5])

t.test(tab2data$naturaltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$naturaltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$naturaltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$biologisttr)
sd(tab2data$biologisttr)
mean(tab2data$biologisttr[tab2data$pidrep<3.5])
sd(tab2data$biologisttr[tab2data$pidrep<3.5])
mean(tab2data$biologisttr[tab2data$pidrep>4.5])
sd(tab2data$biologisttr[tab2data$pidrep>4.5])

t.test(tab2data$biologisttr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$biologisttr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$biologisttr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$chemisttr)
sd(tab2data$chemisttr)
mean(tab2data$chemisttr[tab2data$pidrep<3.5])
sd(tab2data$chemisttr[tab2data$pidrep<3.5])
mean(tab2data$chemisttr[tab2data$pidrep>4.5])
sd(tab2data$chemisttr[tab2data$pidrep>4.5])

t.test(tab2data$chemisttr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$chemisttr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$chemisttr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$physictr)
sd(tab2data$physictr)
mean(tab2data$physictr[tab2data$pidrep<3.5])
sd(tab2data$physictr[tab2data$pidrep<3.5])
mean(tab2data$physictr[tab2data$pidrep>4.5])
sd(tab2data$physictr[tab2data$pidrep>4.5])

t.test(tab2data$physictr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$physictr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$physictr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$socscientiststr <- rowMeans(subset(tab2data, select=c(psychtr, econtr, poliscitr)), na.rm = TRUE)

mean(tab2data$socscientiststr)
sd(tab2data$socscientiststr)
mean(tab2data$socscientiststr[tab2data$pidrep<3.5])
sd(tab2data$socscientiststr[tab2data$pidrep<3.5])
mean(tab2data$socscientiststr[tab2data$pidrep>4.5])
sd(tab2data$socscientiststr[tab2data$pidrep>4.5])

t.test(tab2data$socscientiststr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$socscientiststr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$socscientiststr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$psychtr)
sd(tab2data$psychtr)
mean(tab2data$psychtr[tab2data$pidrep<3.5])
sd(tab2data$psychtr[tab2data$pidrep<3.5])
mean(tab2data$psychtr[tab2data$pidrep>4.5])
sd(tab2data$psychtr[tab2data$pidrep>4.5])

t.test(tab2data$psychtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$psychtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$psychtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$econtr)
sd(tab2data$econtr)
mean(tab2data$econtr[tab2data$pidrep<3.5])
sd(tab2data$econtr[tab2data$pidrep<3.5])
mean(tab2data$econtr[tab2data$pidrep>4.5])
sd(tab2data$econtr[tab2data$pidrep>4.5])

t.test(tab2data$econtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$econtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$econtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$poliscitr)
sd(tab2data$poliscitr)
mean(tab2data$poliscitr[tab2data$pidrep<3.5])
sd(tab2data$poliscitr[tab2data$pidrep<3.5])
mean(tab2data$poliscitr[tab2data$pidrep>4.5])
sd(tab2data$poliscitr[tab2data$pidrep>4.5])

t.test(tab2data$poliscitr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$poliscitr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$poliscitr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$climatetr)
sd(tab2data$climatetr)
mean(tab2data$climatetr[tab2data$pidrep<3.5])
sd(tab2data$climatetr[tab2data$pidrep<3.5])
mean(tab2data$climatetr[tab2data$pidrep>4.5])
sd(tab2data$climatetr[tab2data$pidrep>4.5])

t.test(tab2data$climatetr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$climatetr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$climatetr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$covidtr)
sd(tab2data$covidtr)
mean(tab2data$covidtr[tab2data$pidrep<3.5])
sd(tab2data$covidtr[tab2data$pidrep<3.5])
mean(tab2data$covidtr[tab2data$pidrep>4.5])
sd(tab2data$covidtr[tab2data$pidrep>4.5])

t.test(tab2data$covidtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$covidtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$covidtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

# Perception Scores
mean(tab2data$scimethodim)
sd(tab2data$scimethodim)
mean(tab2data$scimethodim[tab2data$pidrep<3.5])
sd(tab2data$scimethodim[tab2data$pidrep<3.5])
mean(tab2data$scimethodim[tab2data$pidrep>4.5])
sd(tab2data$scimethodim[tab2data$pidrep>4.5])

mean(tab2data$biologim)
sd(tab2data$biologim)
mean(tab2data$biologim[tab2data$pidrep<3.5])
sd(tab2data$biologim[tab2data$pidrep<3.5])
mean(tab2data$biologim[tab2data$pidrep>4.5])
sd(tab2data$biologim[tab2data$pidrep>4.5])

t.test(tab2data$biologim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$biologim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$biologim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$chemim)
sd(tab2data$chemim)
mean(tab2data$chemim[tab2data$pidrep<3.5])
sd(tab2data$chemim[tab2data$pidrep<3.5])
mean(tab2data$chemim[tab2data$pidrep>4.5])
sd(tab2data$chemim[tab2data$pidrep>4.5])

t.test(tab2data$chemim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$chemim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$chemim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$physicistim)
sd(tab2data$physicistim)
mean(tab2data$physicistim[tab2data$pidrep<3.5])
sd(tab2data$physicistim[tab2data$pidrep<3.5])
mean(tab2data$physicistim[tab2data$pidrep>4.5])
sd(tab2data$physicistim[tab2data$pidrep>4.5])

t.test(tab2data$physicistim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$physicistim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$physicistim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")
